Contents

1 Introduction

This tutorial walks through the same dataset used in the “multiome tutorial - archR workflow”. This is a dataset generated by infecting LNCaP cells with NKX2-1 and GATA6 to examine the effects of these TFs on AR activity.

2 Installation

Epiregulon is currently available on R/dev

library(epiregulon)

Alternatively, you could install from gitlab

devtools::install_github(repo ='xiaosaiyao/epiregulon')

library(epiregulon)

3 Data preparation

Single cell preprocessing needs to performed by user’s favorite methods prior to using Epiregulon. The following components are required:
1. Peak matrix from scATAC-seq
2. Gene expression matrix from either paired or unpaired scRNA-seq. RNA-seq integration needs to be performed for unpaired dataset.
3. Dimensionality reduction matrix from with either single modalities or joint scRNA-seq and scATAC-seq

Multiome data can now be conveniently processed by initiate.archr and then gp.sa.archr to obtain peak matrices. Finally, the archR project can be uploaded into DatasetDB as a MultiAssayExperiment object using maw.archr::importArchr or maw.archr::create.mae.with.multiple.sces.from.archr

# load the MAE object
library(SingleCellExperiment)
mae <- dsassembly::getDataset("DS000013080")
#> 'version=' not specified, using the latest version (2) instead
#> Error in read_csv(path, is_compressed = identical(compression, "gzip"),  : 
#>   encountered empty line in a file with non-zero columns
#> falling back to 'read.csv' for a malformed CSV data frame

# peak matrix
PeakMatrix <- mae[["PeakMatrix"]]

# expression matrix
GeneExpressionMatrix <- mae[["GeneExpressionMatrix"]]
rownames(GeneExpressionMatrix) <- rowData(GeneExpressionMatrix)$name

# dimensional reduction matrix
reducedDimMatrix <- reducedDim(mae[['TileMatrix500']], "LSI_ATAC")

Visualize singleCellExperiment by UMAP

# transfer UMAP_combined from TileMatrix to GeneExpressionMatrix
reducedDim(GeneExpressionMatrix, "UMAP_Combined") <- reducedDim(mae[['TileMatrix500']], "UMAP_Combined")
scater::plotReducedDim(GeneExpressionMatrix, 
                       dimred = "UMAP_Combined", 
                       text_by = "Clusters", 
                       colour_by = "Clusters")

4 Quick start

4.1 Retrieve bulk TF ChIP-seq binding sites

First, we retrieve the information of TF binding sites collected from Cistrome and ENCODE ChIP-seq, which are hosted on Genomitory. Currently, human genomes HG19 and HG38 and mouse mm10 are available.

grl <- getTFMotifInfo(genome = "hg38")
#> redirecting from 'GMTY162:hg38_motif_bed_granges@REVISION-4' to 'GMTY162:hg38_motif_bed_granges@24c22e4f4204fd48f2133374e6c1abf086a017d5'
head(grl)
#> GRangesList object of length 6:
#> $`5-hmC`
#> GRanges object with 24048 ranges and 0 metadata columns:
#>           seqnames            ranges strand
#>              <Rle>         <IRanges>  <Rle>
#>       [1]     chr1       10000-10685      *
#>       [2]     chr1       13362-13694      *
#>       [3]     chr1       29631-29989      *
#>       [4]     chr1       40454-40754      *
#>       [5]     chr1     135395-135871      *
#>       ...      ...               ...    ...
#>   [24044]     chrY 56864377-56864627      *
#>   [24045]     chrY 56876124-56876182      *
#>   [24046]     chrM           84-2450      *
#>   [24047]     chrM       13613-14955      *
#>   [24048]     chrM       15134-16490      *
#>   -------
#>   seqinfo: 25 sequences from an unspecified genome; no seqlengths
#> 
#> ...
#> <5 more elements>

4.3 Add TF motif binding to peaks

The next step is to add the TF binding information by overlapping regions of the peak matrix with the bulk chip-seq database loaded in 2. The user can supply either an archR project path and this function will retrieve the peak matrix, or a peakMatrix in the form of a Granges object or RangedSummarizedExperiment.


overlap <- addTFMotifInfo(grl = grl, p2g = p2g, peakMatrix = PeakMatrix)
#> Computing overlap...
#> Success!
head(overlap)
#>   idxATAC idxTF       tf
#> 1       1     2     5-mC
#> 2       1    22 AML1-ETO
#> 3       1    25       AR
#> 4       1    49     ATF1
#> 5       1    50     ATF2
#> 6       1    51     ATF3

4.4 Generate regulons

A long format dataframe, representing the inferred regulons, is then generated. The dataframe consists of three columns:

  • tf (transcription factor)
  • target gene
  • peak to gene correlation between tf and target gene

regulon <- getRegulon(p2g = p2g, overlap = overlap, aggregate = FALSE)
head(regulon)
#>   idxATAC idxTF       tf  chr  start    end idxRNA    target      corr distance       FDR
#> 1       1     2     5-mC chr1 817121 817621     15 LINC01409 0.5368078    38164 0.3395926
#> 2       1    22 AML1-ETO chr1 817121 817621     15 LINC01409 0.5368078    38164 0.3395926
#> 3       1    25       AR chr1 817121 817621     15 LINC01409 0.5368078    38164 0.3395926
#> 4       1    49     ATF1 chr1 817121 817621     15 LINC01409 0.5368078    38164 0.3395926
#> 5       1    50     ATF2 chr1 817121 817621     15 LINC01409 0.5368078    38164 0.3395926
#> 6       1    51     ATF3 chr1 817121 817621     15 LINC01409 0.5368078    38164 0.3395926

4.6 Add Weights

While the `pruneRegulon’ function provides statistics on the joint occurrence of TF-RE-TG, we would like to further estimate the strength of regulation. Biologically, this can be interpreted as the magnitude of gene expression changes induced by transcription factor activity. Epiregulon estimates the regulatory potential using one of the four measures: 1) correlation between TF and target gene expression, 2) mutual information between the TF and target gene expression, 3) Wilcoxon test statistics of target gene expression in cells jointly expressing all 3 elements vs cells that do not, or 4) log 2 fold difference of target gene expression in cells jointly expressing all 3 elements vs cells that do not.

Three measures (correlation, Wilcoxon statistics and log 2 fold difference) give both the magnitude and directionality of changes whereas mutational information is always positive. The correlation and mutual information statistics are computed on the grouped pseudobulks by user-supplied cluster labels, whereas the Wilcoxon and log fold change group cells based on the joint expression of TF, RE and TG in each single cell.


regulon.w <- addWeights(regulon = pruned.regulon,
                        expMatrix  = GeneExpressionMatrix,
                        exp_assay  = "counts",
                        peakMatrix = PeakMatrix,
                        peak_assay = "counts",
                        clusters = GeneExpressionMatrix$Clusters,
                        block_factor = NULL,
                        tf_re.merge = TRUE,
                        method = "corr")
head(regulon.w)
#>       idxATAC idxTF   tf  chr   start     end idxRNA     target      corr distance
#> 7916       67     5 ADNP chr1 1207926 1208426     52      ACAP3 0.5467671    92832
#> 12142     121     5 ADNP chr1 1375168 1375668     56       CPTP 0.5023777    50167
#> 25229     502     5 ADNP chr1 6244982 6245482    173     RNF207 0.6336901    38664
#> 28480     631     5 ADNP chr1 7961380 7961880    210 AL034417.4 0.5311411    27252
#> 30872     655     5 ADNP chr1 8061169 8061669    208    TNFRSF9 0.8219272   120104
#> 30873     655     5 ADNP chr1 8061169 8061669    210 AL034417.4 0.5263050    69836
#>               FDR     pval_all    pval_C1 pval_C2   pval_C3      pval_C4   pval_C5
#> 7916  0.316224023 0.0322720607 0.03803691     NaN       NaN 7.999728e-01 0.4893080
#> 12142 0.421315952 0.0376819525 0.25965219     NaN 0.2116851 2.050637e-01 0.1152127
#> 25229 0.137379482 0.7939385141 0.02686561     NaN 0.7218053 3.715971e-01 0.7645890
#> 28480 0.353112324 0.1125311036 0.41544389     NaN       NaN 6.932539e-01 0.9722747
#> 30872 0.002547118 0.0005790391 0.10494073     NaN       NaN 6.168540e-19 0.8212453
#> 30873 0.365152209 0.1376811319 0.77532454     NaN       NaN 8.751367e-01 0.7299816
#>            pval_C6   stats_all   stats_C1 stats_C2  stats_C3    stats_C4    stats_C5
#> 7916  0.0020815333  4.58398600 4.30337079      NaN       NaN  0.06420262 0.478047626
#> 12142 0.0528806959  4.31932114 1.27061046      NaN 1.5598641  1.60593817 2.481222730
#> 25229 0.2570916475  0.06822528 4.89942722      NaN 0.1267695  0.79831689 0.089676389
#> 28480 0.0396182359  2.51831034 0.66316759      NaN       NaN  0.15558543 0.001207944
#> 30872 0.8848118526 11.84218276 2.62880407      NaN       NaN 79.01354997 0.051050531
#> 30873 0.0001881518  2.20368640 0.08146298      NaN       NaN  0.02469211 0.119128520
#>          stats_C6 padj_all padj_C1 padj_C2 padj_C3      padj_C4 padj_C5 padj_C6    weight
#> 7916   9.47621193        1       1     NaN     NaN 1.000000e+00       1       1 0.8723337
#> 12142  3.74768523        1       1     NaN       1 1.000000e+00       1       1 0.6461784
#> 25229  1.28435077        1       1     NaN       1 1.000000e+00       1       1 0.7622672
#> 28480  4.23415947        1       1     NaN     NaN 1.000000e+00       1       1 0.9617845
#> 30872  0.02098793        1       1     NaN     NaN 2.163537e-12       1       1 0.8466490
#> 30873 13.94585792        1       1     NaN     NaN 1.000000e+00       1       1 0.4693631

4.7 Calculate TF activity

Finally, the activities for a specific TF in each cell are computed by averaging expressions of target genes linked to the TF weighted by the test statistics of choice, chosen from either correlation, mutual information, Wilcoxon test statistics or log fold change. \[y=\frac{1}{n}\sum_{i=1}^{n} x_i * weights_i\] where \(y\) is the activity of a TF for a cell \(n\) is the total number of targets for a TF \(x_i\) is the log count expression of target i where i in {1,2,…,n} \(weights_i\) is the weight of TF and target i

score.combine <- calculateActivity(expMatrix = GeneExpressionMatrix, 
                                   regulon = regulon.w, 
                                   mode = "weight", 
                                   method = "weightedMean", 
                                   exp_assay = "counts",
                                   normalize = FALSE)
#> calculating TF activity from regulon using weightedmean

4.8 Perform differential activity


markers <- findDifferentialActivity(activity_matrix = score.combine, 
                                    groups = GeneExpressionMatrix$Clusters, 
                                    pval.type = "some", 
                                    direction = "up", 
                                    test.type = "t")

Take the top TFs

markers.sig <- getSigGenes(markers, topgenes = 5 )
#> Using a logFC cutoff of 0.3 for class C1
#> Using a logFC cutoff of 0.1 for class C2
#> Using a logFC cutoff of 0.2 for class C3
#> Using a logFC cutoff of 0.1 for class C4
#> Using a logFC cutoff of 0.2 for class C5
#> Using a logFC cutoff of 0 for class C6

4.9 Visualize the results

First visualize the known differential TFs by bubble plot

plotBubble(activity_matrix = score.combine, 
           tf = c("NKX2-1","GATA6","FOXA1","FOXA2", "AR"), 
           clusters = GeneExpressionMatrix$Clusters)

Then visualize the most differential TFs by clusters

plotBubble(activity_matrix = score.combine, 
           tf = markers.sig$tf, 
           clusters = GeneExpressionMatrix$Clusters)

Visualize the known differential TFs by violin plot. Note there is no activity calculated for SOX2 because the expression of SOX2 is 0 in all cells.

plotActivityViolin(activity_matrix = score.combine, 
                   tf = c("NKX2-1","GATA6","FOXA1","FOXA2", "AR", "SOX2"), 
                   clusters = GeneExpressionMatrix$Clusters)
#> SOX2 not found in activity matrix. Excluded from plots

Visualize the known differential TFs by UMAP


plotActivityDim(sce = GeneExpressionMatrix, 
                activity_matrix = score.combine, 
                tf = c("NKX2-1","GATA6","FOXA1","FOXA2", "AR", "SOX2"), 
                dimtype = "UMAP_Combined", 
                label = "Clusters", 
                point_size = 1, 
                ncol = 3)
#> SOX2 not found in activity matrix. Excluded from plots

Visualize the gene expression of the regulons by heatmap

plotHeatmapRegulon(sce=GeneExpressionMatrix,
                   tfs=c("GATA6","NKX2-1"),
                   regulon=regulon.w,
                   regulon_cutoff=0.1,
                   downsample=1000,
                   cell_attributes="Clusters",
                   col_gap="Clusters",
                   exprs_values="counts",
                   name="regulon heatmap")

plotHeatmapActivity(activity=score.combine,
                    sce=GeneExpressionMatrix,
                    tfs=rownames(score.combine),
                    downsample=5000,
                    cell_attributes="Clusters",
                    col_gap="Clusters",
                    name = "Activity")

## Geneset enrichment

Sometimes we are interested to know what pathways are enriched in the regulon of a particular TF. We can perform geneset enrichment using the enricher function from clusterProfiler.

#retrieve genesets
H <- EnrichmentBrowser::getGenesets(org = "hsa", db = "msigdb", cat = "H", 
                                    gene.id.type = "SYMBOL" )
C6 <- EnrichmentBrowser::getGenesets(org = "hsa", db = "msigdb", cat = "C6",
                                     gene.id.type = "SYMBOL" )

#combine genesets and convert genesets to be compatible with enricher
gs <- c(H,C6)
gs.list <- do.call(rbind,lapply(names(gs), function(x) 
  {data.frame(gs=x, genes=gs[[x]])}))

enrichresults <- regulonEnrich(TF = c("GATA6","AR"), 
                               regulon = regulon.w, 
                               corr = "weight",
                               corr_cutoff = 0.1, 
                               genesets = gs.list)
#> GATA6
#> 
#> AR

#plot results
enrichPlot(results = enrichresults, )

## Network analysis

We can visualize the genesets as a network


plotGseaNetwork(tf = names(enrichresults), 
                enrichresults = enrichresults,
                p.adj_cutoff = 0.1,
                ntop_pathways = 10)

4.10 Differential networks

We are interested in understanding the differential networks between two conditions and determining which transcription factors account for the differences in the topology of networks. The pruned regulons with cluster-specific test statistics computed by pruneRegulon can be used to generate cluster-specific networks based on user-defined cutoffs and to visualize differential networks for transcription factors of interest. In this dataset, the GATA6 gene was only expressed in cluster 1 (C1) and NKX2-1 was only expressed in cluster 3 (C3). If we visualize the target genes of GATA6, we can see that C1 has many more target genes of GATA6 compared to C5, a cluster that does not express GATA6. Similarly, NKX2-1 target genes are confined to C3 which is the only cluster that exogenously expresses NKX2-1.


plotDiffNetwork(pruned.regulon,
                cutoff = 1,
                tf = c("GATA6"),
                groups  = c("stats_C1","stats_C5"),
                layout = "stress")


plotDiffNetwork(pruned.regulon,
                cutoff = 1,
                tf = c("NKX2-1"),
                groups  = c("stats_C3","stats_C5"),
                layout = "stress")

We can also visualize how transcription factors relate to other transcription factors in each cluster.


C1_network <- buildGraph(pruned.regulon[pruned.regulon$stats_C1>1 & pruned.regulon$tf %in% c("GATA6","FOXA1","AR"),], weights = "stats_C1")
C5_network <- buildGraph(pruned.regulon[pruned.regulon$stats_C5>1 & pruned.regulon$tf %in% c("GATA6","FOXA1","AR"),], weights = "stats_C5")

plotEpiregulonNetwork(C1_network,
                        layout = "sugiyama",
                        tfs_to_highlight = c("GATA6","FOXA1","AR")) + ggplot2::ggtitle ("C1")


plotEpiregulonNetwork(C5_network,
                        layout = "sugiyama",
                        tfs_to_highlight = c("GATA6","FOXA1","AR")) + ggplot2::ggtitle ("C5")

To systematically examine the differential network topology between two clusters, we perform an edge subtraction between two graphs, using weights computed by pruneRegulon. We then calculate the degree centrality of the weighted differential graphs and if desired, normalize the differential centrality against the total number of edges. The default normalization function is sqrt as it preserves both the difference in the number of edges (but scaled by sqrt) and the differences in the weights. If the user only wants to examine the differences in the averaged weights, the FUN argument can be changed to identity. Finally, we rank the transcription factors by (normalized) differential centrality.


# rank by differential centrality
C1_network <- buildGraph(pruned.regulon, weights = "stats_C1")
C5_network <- buildGraph(pruned.regulon, weights = "stats_C5")

diff_graph <- buildDiffGraph(C1_network, C5_network)
diff_graph <- addCentrality(diff_graph)
diff_graph <- normalizeCentrality(diff_graph)
rank_table <- rankTfs(diff_graph)

library(ggplot2)
ggplot(rank_table, aes(x = rank, y = centrality)) +
    geom_point() +
    ggrepel::geom_text_repel(data = head(rank_table, 5), aes(label = tf)) +
    theme_classic()

5 Session Info

sessionInfo()
#> R version 4.2.0 (2022-04-22)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.6 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/local/lib/R/lib/libRblas.so
#> LAPACK: /usr/local/lib/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C            LC_COLLATE=C        
#>  [5] LC_MONETARY=C        LC_MESSAGES=C        LC_PAPER=C           LC_NAME=C           
#>  [9] LC_ADDRESS=C         LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] ggplot2_3.4.0               org.Hs.eg.db_3.16.0         AnnotationDbi_1.60.0       
#>  [4] msigdbr_7.5.1               BiocStyle_2.26.0            epiregulon_1.0.22          
#>  [7] SingleCellExperiment_1.20.0 SummarizedExperiment_1.29.1 Biobase_2.58.0             
#> [10] GenomicRanges_1.50.2        GenomeInfoDb_1.34.4         IRanges_2.32.0             
#> [13] S4Vectors_0.36.1            BiocGenerics_0.44.0         MatrixGenerics_1.10.0      
#> [16] matrixStats_0.63.0         
#> 
#> loaded via a namespace (and not attached):
#>   [1] utf8_1.2.2                  tidyselect_1.2.0            RSQLite_2.2.19             
#>   [4] grid_4.2.0                  BiocParallel_1.32.4         scatterpie_0.1.8           
#>   [7] munsell_0.5.0               ScaledMatrix_1.6.0          base64url_1.4              
#>  [10] codetools_0.2-18            statmod_1.4.37              scran_1.26.1               
#>  [13] withr_2.5.0                 GOSemSim_2.24.0             colorspace_2.0-3           
#>  [16] genomitory_2.1.6            filelock_1.0.2              highr_0.9                  
#>  [19] knitr_1.40                  rstudioapi_0.13             DOSE_3.23.3                
#>  [22] labeling_0.4.2              KEGGgraph_1.58.2            GenomeInfoDbData_1.2.9     
#>  [25] polyclip_1.10-4             bit64_4.0.5                 farver_2.1.1               
#>  [28] rhdf5_2.42.0                downloader_0.4              treeio_1.22.0              
#>  [31] vctrs_0.5.1                 generics_0.1.3              gson_0.0.9                 
#>  [34] xfun_0.31                   BiocFileCache_2.6.0         R6_2.5.1                   
#>  [37] doParallel_1.0.17           graphlayouts_0.8.4          ggbeeswarm_0.7.1           
#>  [40] clue_0.3-63                 rsvd_1.0.5                  gp.version_1.5.0           
#>  [43] locfit_1.5-9.6              artificer.matrix_1.3.7      gridGraphics_0.5-1         
#>  [46] fgsea_1.24.0                bitops_1.0-7                rhdf5filters_1.10.0        
#>  [49] cachem_1.0.6                DelayedArray_0.24.0         assertthat_0.2.1           
#>  [52] scales_1.2.1                ggraph_2.1.0                enrichplot_1.18.3          
#>  [55] beeswarm_0.4.0              gtable_0.3.1                beachmat_2.14.0            
#>  [58] Cairo_1.6-0                 metacommons_1.9.0           tidygraph_1.2.2            
#>  [61] rlang_1.0.6                 splines_4.2.0               GlobalOptions_0.1.2        
#>  [64] lazyeval_0.2.2              dsdb.schemas_0.99.1         checkmate_2.1.0            
#>  [67] gp.auth_1.7.0               BiocManager_1.30.19         yaml_2.3.5                 
#>  [70] reshape2_1.4.4              backports_1.4.1             qvalue_2.30.0              
#>  [73] clusterProfiler_4.6.0       tools_4.2.0                 bookdown_0.30              
#>  [76] ggplotify_0.1.0             jquerylib_0.1.4             RColorBrewer_1.1-3         
#>  [79] artificer.base_1.3.19       MultiAssayExperiment_1.24.0 Rcpp_1.0.9                 
#>  [82] plyr_1.8.8                  sparseMatrixStats_1.10.0    zlibbioc_1.44.0            
#>  [85] purrr_0.3.5                 RCurl_1.98-1.9              artificer.ranges_1.3.4     
#>  [88] GetoptLong_1.0.5            viridis_0.6.2               cowplot_1.1.1              
#>  [91] ShadowArray_1.7.1           ggrepel_0.9.2               cluster_2.1.3              
#>  [94] magrittr_2.0.3              data.table_1.14.6           magick_2.7.3               
#>  [97] circlize_0.4.15             patchwork_1.1.2             evaluate_0.19              
#> [100] GSVA_1.46.0                 xtable_1.8-4                HDO.db_0.99.0              
#> [103] XML_3.99-0.13               artificer.schemas_0.99.2    artificer.sce_1.3.4        
#> [106] gridExtra_2.3               shape_1.4.6                 compiler_4.2.0             
#> [109] scater_1.26.1               tibble_3.1.8                shadowtext_0.1.2           
#> [112] crayon_1.5.1                gp.cache_1.7.1              htmltools_0.5.4            
#> [115] ggfun_0.0.9                 aplot_0.1.9                 tidyr_1.2.1                
#> [118] ArtifactDB_1.9.5            DBI_1.1.3                   tweenr_2.0.2               
#> [121] dbplyr_2.2.1                ComplexHeatmap_2.14.0       MASS_7.3-58.1              
#> [124] rappdirs_0.3.3              babelgene_22.9              Matrix_1.5-3               
#> [127] cli_3.4.1                   artificer.mae_1.3.4         genomitory.schemas_0.99.0  
#> [130] parallel_4.2.0              metapod_1.6.0               igraph_1.3.5               
#> [133] pkgconfig_2.0.3             getPass_0.2-2               dsdb.plus_1.3.2            
#> [136] scuttle_1.8.3               foreach_1.5.2               ggtree_3.6.2               
#> [139] annotate_1.76.0             vipor_0.4.5                 bslib_0.4.2                
#> [142] dqrng_0.3.0                 ArchR_1.0.2                 XVector_0.38.0             
#> [145] yulab.utils_0.0.5           stringr_1.4.0               digest_0.6.29              
#> [148] graph_1.76.0                Biostrings_2.66.0           fastmatch_1.1-3            
#> [151] rmarkdown_2.18              tidytree_0.4.1              artificer.se_1.3.5         
#> [154] edgeR_3.40.1                DelayedMatrixStats_1.20.0   GSEABase_1.60.0            
#> [157] curl_4.3.2                  rjson_0.2.21                nlme_3.1-161               
#> [160] lifecycle_1.0.3             jsonlite_1.8.4              Rhdf5lib_1.20.0            
#> [163] dsassembly_1.7.6            BiocNeighbors_1.16.0        viridisLite_0.4.1          
#> [166] limma_3.54.0                fansi_1.0.3                 pillar_1.8.1               
#> [169] lattice_0.20-45             GO.db_3.16.0                KEGGREST_1.38.0            
#> [172] fastmap_1.1.0               httr_1.4.3                  glue_1.6.2                 
#> [175] png_0.1-8                   iterators_1.0.14            bluster_1.8.0              
#> [178] bit_4.0.5                   Rgraphviz_2.42.0            ggforce_0.4.1              
#> [181] stringi_1.7.6               sass_0.4.4                  HDF5Array_1.26.0           
#> [184] BiocBaseUtils_1.1.0         EnrichmentBrowser_2.28.0    blob_1.2.3                 
#> [187] BiocSingular_1.14.0         memoise_2.0.1               dplyr_1.0.10               
#> [190] ape_5.6-2                   irlba_2.3.5.1